Installing and Loading the R package

maptools: to create spatial objects
sf: to import geosptial data
raster: to create raster object
spatstat: to perform spatial point pattern analysis (SPPA)
tmap: to plot thematic maps
tidyverse for data handling

packages = c('maptools', 'sf', 'raster', 'spatstat', 'tmap', 'tidyverse', 'plotly', 'ggthemes')

for (p in packages){
if(!require(p, character.only = T)){
  install.packages(p)
}
  library(p,character.only = T)
}

Importing Geospatial Data

Import shapefile using st_read() of sf package. Output object is in tibble sf object class.

mpsz_sf <- st_read(dsn = 'data/shapefile', layer = 'MP14_SUBZONE_WEB_PL')

Reading layer `MP14_SUBZONE_WEB_PL' from data source 
  `C:\wlwong2018\IS415 Blog\_posts\2021-09-13-in-class-exercise-5\data\shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21

Projection is in SVY21.

Import aspatial data from rds folder

read_rds() of readr package is used instead of readRDS() of base R is used. This is because output of read_rds() is in tibble object.

childcare <- read_rds("data/rds/childcare.rds")
CHAS <- read_rds("data/rds/CHAS.rds")

Note that there are some data issue in childcare data frame because Lat and Lng should be in numeric data type. The coordinate fields seem to be in decimal degree. Hence wgs referencing system is assumed.

Converting the aspatial data frame into sf objects

CHAS_sf <- st_as_sf(CHAS,
                    coords = c("X_COORDINATE", "Y_COORDINATE"),
                    crs = 3414)

Note: st_as_sf() accept coordinates in character data type.
lat long in decimal degrees

childcare_sf <- st_as_sf(childcare,
                    coords = c("Lng", "Lat"),
                    crs = 4326)  %>%
  st_transform(crs = 3414)

Plot to review

alpha: to set transparency

tmap_mode("view")

tm_shape(childcare_sf) +
  tm_dots(alpha = 0.4,
          col = "blue",
          size = 0.05) +
tm_shape(CHAS_sf) +
  tm_dots(alpha = 0.4,
          col = "red",
          size = 0.05)

Geospatial Data Wrangling

Convert sf to Spatial* classes

as_Spatial() of sf package

childcare <- as_Spatial(childcare_sf)
CHAS <- as_Spatial(CHAS_sf)
mpsz <- as_Spatial(mpsz_sf)

Convert Spatial* dataframe into Spatial* objects

as.SpatialPoint() of as.SpatialPolygon() of maptools package

childcare_sp <- as(childcare, "SpatialPoints") # data will be dropped. 
CHAS_sp <- as(CHAS, "SpatialPoints") # data will be dropped. 
mpsz_sp <- as(mpsz, "SpatialPolygons") # data will be dropped.

Convert from Spatial* objects into ppp object

to drop project info. only X, Y coordinates will be left.
as.ppp() of maptools package

childcare_ppp <- as(childcare_sp, "ppp")
CHAS_ppp <- as(CHAS_sp, "ppp")

Remove duplicated points using jitter

to avoid points overlapping

childcare_ppp_jit <- rjitter(childcare_ppp,
                             retry = TRUE,
                             nsim = 1,
                             drop = TRUE)

any(duplicated(childcare_ppp_jit))

[1] FALSE

CHAS_ppp_jit <- rjitter(CHAS_ppp,
                             retry = TRUE,
                             nsim = 1,
                             drop = TRUE)

any(duplicated(CHAS_ppp_jit))

[1] FALSE

Note:

must be sf or sp to plot a tmap.
For now, we can only plot using sf plot()

Extract Punggol Planning Area

Remember to add the comma at the end

pg <- mpsz[mpsz@data$PLN_AREA_N == "PUNGGOL",]

Convert SpatialPolygonsDataFrame into SpatialPolygons object

pg_sp <- as(pg, "SpatialPolygons")

COnvert SpatialPolygons into owin object

owin: an object needed by spatstat to define irregular window

pg_owin <- as(pg_sp, "owin")

Extract spatial points within owin

childcare_pg <- childcare_ppp_jit[pg_owin]
CHAS_pg <- CHAS_ppp_jit[pg_owin]

Plotting with childcare_ppp_jit would not work. Does not show the boundary of Punggol

plot(childcare_pg)

L-function

L_childcare <- envelope(childcare_pg,
                        Lest,
                        nsim=99,
                        rank=1,
                        global = TRUE)

Generating 99 simulations of CSR  ...
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70,
71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.

Done.

Code chunk for plotting interactive L-function

title <- "Pairwise Distance: L function"

Lcsr_df <- as.data.frame(L_childcare)

colour=c("#0D657D","#ee770d","#D3D3D3")
csr_plot <- ggplot(Lcsr_df, aes(r, obs-r))+
  # plot observed value
  geom_line(colour=c("#4d4d4d"))+
  geom_line(aes(r,theo-r), colour="red", linetype = "dashed")+
  # plot simulation envelopes
  geom_ribbon(aes(ymin=lo-r,ymax=hi-r),alpha=0.1, colour=c("#91bfdb")) +
  xlab("Distance r (m)") +
  ylab("L(r)-r") +
  geom_rug(data=Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,], sides="b", colour=colour[1])  +
  geom_rug(data=Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,], sides="b", colour=colour[2]) +
  geom_rug(data=Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,], sides="b", color=colour[3]) +
  theme_tufte()+
  ggtitle(title)

text1<-"Significant clustering"
text2<-"Significant segregation"
text3<-"Not significant clustering/segregation"

# the below conditional statement is required to ensure that the labels (text1/2/3) are assigned to the correct traces
if (nrow(Lcsr_df[Lcsr_df$obs > Lcsr_df$hi,])==0){ 
  if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text3, traces = 4) %>%
      rangeslider() 
  }else if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){ 
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      rangeslider() 
  }else {
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text2, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider() 
  }
} else if (nrow(Lcsr_df[Lcsr_df$obs < Lcsr_df$lo,])==0){
  if (nrow(Lcsr_df[Lcsr_df$obs >= Lcsr_df$lo & Lcsr_df$obs <= Lcsr_df$hi,])==0){
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      rangeslider() 
  } else{
    ggplotly(csr_plot, dynamicTicks=T) %>%
      style(text = text1, traces = 4) %>%
      style(text = text3, traces = 5) %>%
      rangeslider()
  }
} else{
  ggplotly(csr_plot, dynamicTicks=T) %>%
    style(text = text1, traces = 4) %>%
    style(text = text2, traces = 5) %>%
    style(text = text3, traces = 6) %>%
    rangeslider()
  }

In-class Exercise 5

Installing and Loading the R package

Importing Geospatial Data

Import aspatial data from rds folder

Converting the aspatial data frame into sf objects

Plot to review

Geospatial Data Wrangling

Convert sf to Spatial* classes

Convert Spatial* dataframe into Spatial* objects

Convert from Spatial* objects into ppp object

Remove duplicated points using jitter

Extract Punggol Planning Area

Convert SpatialPolygonsDataFrame into SpatialPolygons object

COnvert SpatialPolygons into owin object

Extract spatial points within owin

L-function

Code chunk for plotting interactive L-function